﻿using System;
using System.Collections.Generic;
using System.Text;

namespace VerticalProcesses
{
    public abstract class EvaporationEquations
    {
        public abstract double Evaporate(DateTime TimeStep, double TimeStepHrs, double gruTa, double gruTmin, double gruTmax, double latitude, double longitude, double elevation, double timezone, double slope, double aspect);
    }
    /// <summary>
    /// Calculates Potential Evaporation according to Penman Monteith
    /// This follows the FAO recommendation for using PM when only air temperature 
    /// is available and uses the same constants as Environment Canada. 
    /// IMPORTANT: contrary to FAO and EC, this formulation accounts for slope 
    /// and aspect effects by correcting clearsky solar radiation for slope and aspect
    /// 
    /// </summary>
    public class PenmanMonteith : EvaporationEquations
    {
        public override double Evaporate(DateTime TimeStep, double TimeStepHrs, double gruTa, double gruTmin, double gruTmax, double latitude, double longitude, double elevation, double timezone, double slope, double aspect)
        {
            //Constants
            //These are constants as used by Environment Canada

            double WS = 2;	       //m/s  wind speed 
            double ALB = 0.23;	   //Albedo or canopy reflection coefficient
            double GSC = 0.082;	   //Solar constant
            double KRS = 0.16;     //Interior=0.16; Coastal=0.19	0.16	
            double K0 = 2;	       //Ko = 0 °C for humid and sub-humid climates, Ko = 2 °C for arid and semi-arid climates

            //1. calculate Tmean and Julian Days
            double Tmean = (gruTmin + gruTmax) / 2.0;
            int Jday = TimeStep.DayOfYear;

            //2. Calculate actual vapour pressure, ea [kPa]
            //Use minimum temperature and adjustment factor depending on climate classification humid or semi-arid. 
            double Ea = ClimateEq.CalcSaturationVapourPressure(gruTmin - K0);

            //3.	Calculate saturated vapour pressure for Tmin, e(Tmin) [kPa]
            double Etmin = ClimateEq.CalcSaturationVapourPressure(gruTmin);
            
            //4.	Calculate saturated vapour pressure for Tmax, e(Tmax) [kPa]
            double Etmax = ClimateEq.CalcSaturationVapourPressure(gruTmax);

            
            //5.	Calculate saturated vapour pressure, es [kPa]
            double Es = (Etmin + Etmax) / 2;

            //6.	Calculate inverse relative distance Earth-Sun, dr [rad]
            double dr = RadiationEq.CalcInverseRelativeDistanceEarthSun(Jday);

            //7.	Calculate extraterrestrial radiation, Ra [MJm-2 day-1] 
            double Ra = RadiationEq.CalcExtraterrRad(TimeStep, 24, latitude, longitude, elevation, timezone, slope, aspect);

            //8.	Calculate clear sky solar radiation, Rso [MJm-2 day-1]
            double Rso = RadiationEq.CalcClearSkySolarRadiation(TimeStep, 24, latitude, longitude, elevation, timezone, slope, aspect);

            //9.	Calculate solar radiation, Rs [MJm-2 day-1]
            //Use adjustment factor KRs depending on station location, coastal or interior 
            //where:
            //KRs = 0.16 for interior locations
            //KRs = 0.19 for coastal locations
            double Rs = KRS * Math.Sqrt(gruTmax - gruTmin) * Ra;

            //10.	Calculate net longwave radiation, Rnl [MJm-2 day-1] 
            //where:
            // Bolz   = 4.903 x 10-9 MJK-4m-2day-1
            double Rnl = RadiationEq.CalcDailyNetLongwaveRadiation(gruTmax, gruTmin, Ea, Rs, Rso);

            //11.	Calculate net solar radiation, Rns [MJm-2 day-1]
            double Rns = RadiationEq.CalcNetShortwaveRadiation(0.23, Rs);

            // 12.	Calculate net radiation, Rn [MJm-2 day-1]
            double Rn = RadiationEq.CalcNetRadiation(Rns, Rnl);

            // 13.	Calculate slope vapour pressure,   [kPa  C-1]
            double SLVP = ClimateEq.CalcSlopeOfSaturationPressureCurve(Tmean);

            //14.	Calculate atmospheric pressure, PZ [kPa]
            double PZ = ClimateEq.CalcAtmosphericPressure(elevation, Tmean);

            // 15.	Calculate psychometric constant,   [kPa  C-1]
            double Gamma = ClimateEq.CalcPsyConst(PZ);

            //16.	Calculate evapotranspiration,  ETr
            double ETr = (0.408 * SLVP * Rn + Gamma * 900 / (Tmean + Constants.KelvinScaleOffset) * WS * (Es - Ea)) / (SLVP + Gamma * (1 + 0.34 * WS));
            return Math.Max(ETr, 0.0);
        }
    }
    public class Hamon : EvaporationEquations
    {

        /// <summary>
        /// Calculates potential evaporation according to Hamon (1961)
        /// retuns 0.0 during nighttime
        /// </summary>
        /// <param name="TimeStep"></param>
        /// <param name="gruTa"></param>
        /// <returns>potential evaporation</returns>
        public override double Evaporate(DateTime TimeStep, double TimeStepHrs, double gruTa, double gruTmin, double gruTmax, double latitude, double longitude, double elevation, double timezone, double slope, double aspect)
        {

            //if (gruShortRad == 0.0) // no evaporation during nighttime
            //    return 0.0;
            //// scale evaporation with temperature by calculating ETHamon 
            // as if for a day and prorate relative to daylength
            double P = Math.Asin(.39795 * Math.Cos(.2163108 + 2 * Math.Atan(.9671396 * Math.Tan(.00860 * (TimeStep.DayOfYear - 186)))));
            double dayLength = 24.0 - (24.0 / Math.PI) * Math.Acos((Math.Sin(0.8333 * Math.PI / 180.0)
                + Math.Sin(latitude * Math.PI / 180.0) * Math.Sin(P)) / (Math.Cos(latitude * Math.PI / 180.0) * Math.Cos(P)));  //Approximate Day length
            //double SVD1 = 4.95 * Math.Exp(0.062 * gruTa);-->joels equation
            double SVD = (1324.0 / (gruTa + 273.2) * Math.Exp(17.27 * gruTa / (gruTa + 237.3)));//Dingman D-7, D8c
            return 0.55 * Math.Pow(dayLength / 12.0, 2) * SVD / 100.0 * 25.4 * TimeStepHrs / 24.0;
        }

    }
   
}
